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Abstract. On the basis of a recently proposed strategy of finite element 
integration in time domain for partial differential equations with a singular 
source term, we present a fourth order algorithm for non-rotating black hole 
perturbations in the Regge- Wheeler gauge. Herein, we address even perturbations 
induced by a particle plunging in. The forward time value at the upper node of 
the (r* , t) grid cell is obtained by an algebraic sum of i) the preceding node 
values of the same cell, ii) analytic expressions, related to the jump conditions 
on the wave function and its derivatives, iii) the values of the wave function at 
adjacent cells. In this approach, the numerical integration does not deal with the 
source and potential terms directly, for cells crossed by the particle world line. 
This scheme has also been applied to circular and eccentric orbits and it will be 
object of a forthcoming publication. 



PACS numbers: 04.25.Nx, 04.30.Db, 04.30.Nk, 04.70.Bw, 95.30.Sf 



1. Introduction 

In the scenario of the capture of compact objects by a supermassive black hole of mass 
M, the seized object is compared to a small mass m (henceforth the particle or the 
source) perturbing the background spacetime curvature and generating gravitational 
radiation. A comprehensive introduction to the general relativistic issues related to 
EMRI (Extreme Mass Ratio Inspiral) sources is contained in a topical volume [l]. 

Schwarzschild-Droste (SD) (2j|5] (see Rothman 6| for a justification of this 
terminology), black hole perturbations have been hugely developed in the Regge- 
Wheeler (RW) gauge, before in vacuum [t] and after in the presence of a particle by 
Zerilli [8}fTT] . The first finite difference scheme in time domain was proposed by Lousto 



and Price [12 . The initial conditions, reflecting the past motion of the particle and the 



initial amount of gravitational waves, were parametrised by Martel and Poisson 13 . 

If the gravitational radiation emitted and the mass of the captured object are to 
be taken into account for the determination of the motion of the latter, it is necessary 
to compute the derivatives of the perturbations that implies the third derivative of the 



wave function \l/(r*, i), see e.g. 14 . For a given accuracy 0{h) of the third derivative 
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of the error on "if itself should be ©(/i^). Effectively, the reminder ought to be 0{h^) 
due to the presence in the mesh of the particle that lowers by one more degree the 
convergence order of the code for geometrical effects 15 . We have therefore developed 
a fourth order scheme. 

The complexity in assessing the continuity of the perturbations at the position 
of the particle and the compatibility of the self-force to the harmonic (Lorenz-de 
Dondei[|]) gauge 16 17 has led researchers to convey their efforts to this gauge, as 
commenced by Barack and Lousto 19 1. Conversely, work in harmonic gauge is made 
cumbersome by the presence of a system of ten coupled equations which replace the 
single wave equation of the RW gauge. 

We have proposed 20 21 a finite element method of integration , in RW gauge. 



based on the jump conditions that the wave function and its derivatives have to satisfy 
for the SD black hole perturbations to be continuous at the position of the particle. 
We first deal with the radial trajectory and the associated even parity perturbations, 
while in a forthcoming paper we shall present the circular and eccentric orbital cases, 
referring thus to both odd and even parity perturbations. 

The main feature of this method consists in avoiding the direct and explicit 
integration of the wave equation (the potential and the source term with the 
associated singularities) whenever the grid cells are crossed by the particle. Indeed, 
the information on the wave equation is implicitly given by the jump conditions on 
the wave function and its derivatives. Conversely, for cells not crossed by the particle 
world line, the integrating method might retain the previous approach by Lousto 22 
and Haas 15 . Among the efforts using jump discontinuities, although in a different 



context, it is worthwhile to mention those of Haas 15 , Sopuerta and coworkers 23 - 25 



getting the self-force in a scalar case. For the geodesic gravitational case, like Sopuerta 
and coworkers, 



Zumbusch 
Evans [30 



28 



Jung et al. 
Field et al. 



26 



29 



Chakraborty et al. 27 



rely on spectral methods; 
use a discontinuous Galerkin method; Hopper and 
work partially in frequency domain. Among recent results not based 
on jump discontinuities but concerning fourth order time domain codes, the one 
proposed by Thornburg 31 deals with and adaptive mesh refinement, while Nagar 



and coworkers replace the delta distribution with a narrow Gaussian 32 33 



For the computation of the back-action, this method ensures a well behaved wave 
function at the particle position, since the approach is governed by the analytical 
values of the jump conditions at the particle position. 

In 21 we have provided waveforms at infinity and the wave function at the 
position of the particle at first order. Herein, we focus instead on the improvement 



of the algorithm at fourth order and refer to 21 for all complementary information 



The features of this method can be summarised as follows: 

• Avoidance of direct and explicit integration of the wave equation (the potential 
and the source term with the associated singularities) for the grid cells crossed 
by the particle. 

• Improvement of the reliability, since analytic expressions partly replace numerical 



ones (the replacement is total at first order 20 21 



• Applicability of the method to generic orbits, assuming that the even and odd 
wave equations are satisfied by 5*, respectively R, being 

§ FitzGerald is considered to have also identified the harmonic gauge [18| . 

II A continuity class element, like a Heaviside step distribution, may be seen as an element 

which after integration transforms into an element belonging to the C" class of functions. 



Fourth order integration method 



3 



Geometric units (G — c — 1) are used, unless stated otherwise. The metric 
signature is (— , +, +, +). 



2. The wave equation 



The wave function (its dimension is such that the energy is proportional to '^^ dt), 
in the Moncrief form ^34j and RW gauge Jj , is defined by 



A + 1 



2M 



Xr + 3M 



dr 



(1) 



where K{t,r) and H2{t,r) are the perturbations, and the Zerilli |9 normalisation is 
used for 4'/ . The wave equation is given by the operator Z acting on the wave function 

Z^\t,r) = dl,¥{t,r)-d^t-^\t,r)-V\r)¥{t,r) = S\t,r) , (2) 

where r* ~ r + 2M\n{r/2M — 1) is the tortoise coordinate and the potential V^{r) is 



y'(r)= 1- 



2M \ 2A2 (A + l)r3+6A2Mr2 + ISAM^r + ISM^ 



r- 



r3(Ar+3M)2 



(3) 



being A — 1/2{1 — + 2). The source S\t,r) includes the derivative of the Dirac 
distribution (the latter appear in the process of forming the wave equation out of the 
ten linearised Einstein equations) 



5' 



r(7^-2M) 



2{r - 2M)k 
~ r2(A + l)(Ar + 3M) 

ru{t)] 



r{X + 1) - 3Af 
21P 



3MU°{r ~ 2Mf 
r{Xr + 3A/) 



5[r^r^{t)] \ , (4) 



= E/{\ — 2M/r„) being the time component of the 4- velocity, E = \Jl — 2M/r„o 
the conserved energy per unit mass, and k = AruyJ (21 + l)7r. The geodesic in the 
unperturbed SD metric 2:„(t) = {t„(T), r„(r), 0„(t), (/)„(t)} assumes different forms 
according to the initial conditions. For radial infall of a particle starting from rest at 
finite distance ruo-, fu(t) is the inverse function in coordinate time t of the trajectory 
in the background field, corresponding to the geodesic in proper time r (u stands for 
unperturbed) 



2M 



ruO 



"dJ \2M 



ruo \2MJ \2MJ 



2arctanh 



2M 

r„ 



2M\ 

ruO 



1 



2M 



2M 

ruO 



4M 
ruO 



V2M/ 



3/2 



arctan 



ruQ 



- 1 



(5) 



The expressions above correspond to those in 14 , where some of the errors of 
previously published literature on radial fall are indicated. 
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3. Jump conditions 

From the visual inspection of the Zerilli wave equation ([2]), it is evinced that the 
wave function ^ is of continuity class (the second derivative of the wave function 
is proportional to the first derivative of the Dirac distribution, in itself a class 
element). Thus, the wave function and its derivatives can be written as (the I index 
is dropped henceforth for simplicity of notation) 



(6) 
(7) 
(8) 
(9) 

^-)S',{10) 
)S' , (11) 

where in shortened notation Qi = <d[i — ru{t)], and Q2 = 'd[ru{t)—r] are two 
Heaviside step distributions, while S — 6[r — ru{t)] and S' = — r„(i)] are the 
Dirac delta - and its derivative - distributions. The dot and the prime indicate time 
and space derivatives, respectively. 





= *+ei4 


*-e2 , 










-*;e2 + 


(*+ - 5 , 








-*^e2- 










+ *;;re2 


+ 2 (*+ - s + {^+ - *-) S' , 








^ Ktt(^2 - 


- 2r„ {^+ - *^) (5 - f„ (*+ -^-)S + 


rl (*+ - 






f *:i,e2 


+ (*+ -^:t)S- r„ {^+ - *;) S - f„ 


(*+ - 



3.1. Jump conditions from the wave equation 

For the computation of back-action effects, we need first order derivatives of the 
perturbations and thus third order wave function derivatives. To this end, we operate 
directly on the wave equation, Eq. [2] The source term is cast in the following form 

S{t, r) = G{t, r)5 + F{t, r)6' = G,„(t)5 + ^^,„(t)<5' , (12) 

where G^^ (t) = G^^^ {t)~ -^l (t) '^^^ '^^ properties of the Dirac delta distribution, 
namely 4){r)5' [r — ru{i)\ = (f>r^(^t)^' [r — fuit)] ~ 0r„(t)^ V ~ 1 been used at 

the position of the particle. The determination of the jump conditions imposes the 
transformation of Eq. [2] into the corresponding equation in (r,t) domain (the tortoise 
coordinate can't be inverted). Turning to the r variable, we get (/ = 1 — 2AI/r) 



dl,^ = ffdr<i> + fdl^ 



[//'*+ + fH+„] Gi + + fH-^] 82 + //' (*^ 

+ 2/2 (*+ - *-) (5 + f (*+ - *-) 5' , 



(13) 



+ rl (*+-*-) 5' , 



(14) 



V* = i/^+ei + v*-e2 . (15) 

The notation [^P] stands for the difference {^~^ — ^~)^ and a likewise notation 
is used for the derivatives at the point r„. Equating the coefficients of S', and owing 
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to the above mentioned property of the delta derivative for which (^'+ — '^~)6' = 
[9] 6' — [9^r] we get the jump condition for ^ 

m = 7^^- • (16) 

Equating the coefficients of S, we get the jump condition on the space derivative 

1 



+ (/r„/;„ - ru) m - 2ru-f- m 



(17) 



and therefore the jump condition on the first time derivative 

[*t]=r„^[*]-f„[*.] . (18) 
Since Z^"^ = 0, the coefficients of 6i and ©2 ought to be equal. We thus obtain 

- /r„ - fl [^,rr] + [*] = , (19) 

which is an equation with two unknowns. We circumvent the difficulty by using 
i) the commutativity of the derivatives, [^,tr] = [^,rt]) ii) the transformation 
d/dt = fud/dru, and write 

dru dr^^ 

The jump condition on the second space derivative can now be expressed by 

= J2^^ {^u^^ -rl^^ [*,^] - frJL [*.^] + [*]} • (21) 

The other second derivatives are obtained by 

[*,tr] = = ^ [*,r] - r„ [^,rr] , (22) 

M = |[*,t]-r„ . (23) 

For the third order derivatives, we derive the wave equation with respect to r and 
obtain 



= r.,. ^ -fl— + rl [^,rr] . (20) 



rfj 



rl- f?^\ dru dru' 

+ {frl + frX: - + ^frJi ,rr] " K„ [*] | , (24) 

while deriving with respect to t, we obtain 
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= [*,trt] = [*,rtt] = ^ [*,«] - , (26) 

[*,trr] = [*,rtr] = [*,rrt] = [*,tr] - r"^ [*,«r] , (27) 

urn 

[*,rrr] = [*,,.r] " . (28) 

Finally, we similarly proceed for the fourth derivatives 

U J Vu 

[*,tttr] = [*,ttrt] = [*,frtt] = [^-.rttt] ^ [*,ttt] - [*,t«t] , (30) 

= ^M-r-M*t«.] , (31) 

[*,frrr] = [*,rtrr] = [*,rrtr] [^-.rrrt] [*,trr] " [*,ftrr] , (32) 

clru 

[*,rrrr] = [*,rTr] - V'^ [*,rrrt] • (33) 

uru 

3.1.1. Jump conditions in explicit form. We list hereafter the jump conditions in 
explicit form. 



Jump conditions 



m = ^ (34) 

^ ^ (A + 1)(3M + Ar„) ^ ' 



First derivative jump conditions 

, ^ KEruT-u , 

^ (2M-r„)(3M + Ar„) ^ ' 

^ nE [6M^ + 3A^Ar„ + A(A + l)rg] 
^ ■'■J (A + l)(2M-r„)(3M + Ar„)2 ^""^^ 

Second derivative jump conditions 

[3M3(5A - 3) + 6M2A(A - 3)r„ + 3A^A2(A - l)r2 - 2A2(A + l)r3] 

^^'"^ ~ (A + l)(2M-r„)2(3M + Ar„)3 ^^^^ 

Kig(3M^+3MAr.-Arg)r. 
- (2M-r„)2(3M + Ar.)2 ^^^^ 

''^■"'°- .„M3M + ..A) 
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Third derivative jump conditions 

kE 



r„ (A + 1) (2M - r„)^(3M + r„A)* 



81 (A + 1) + 9r„ (19A2 + 18E'^X+ 

3A + 18^;^) + 9rlX (TA^ + 24S2a - 14A + 2'iE'^ + 3) + 3r^A^ (TA^ 
+36£;2a - llA + 36£2 + 18) + 3r^A3 {SE'^X - 7X + SE"^ - l) M + 

2rlX^ (A + 1) (^^A + 3) (40) 



-nEfu 



r„(2M-r„)^(3M + r„A)^ 



27M'' + 6r„ (5A + 9E=^ - 3) + 3r^A (5A+ 



18i;2 _ 6) + er^A^ (3i;2 - 2) M + 2r^A^ (E^A + l) 



(41) 



kE 



Fourth 

\^ ^rrrr\ 



r3(2M-r„) (3M + r„A) 
12^;^ - 13) M + 2r^A2 [E'^ - l) 



39M^ + 9r„ (3A + 2^;^ - 2) + A (4A+ 



r3 (2M-r„) (3M + r„A) 

derivative jump conditions 
-3kE 



9M2+2r„ (2A + 3E^ - 2) M+2r^A (E^ - l) 



(42) 



(43) 



r2(A + l)(2M-r„)*(3M + r„A)'' 



567 (A + 1) + 162r„ (A + 1) (6A 



+16E^ -5)M^ + 6rl (l39A^ + 738i;^A^ - 123A^ + 162£;^A + 441£;2a 
-171A + 162^;* - 297E'^ + 27) + 12r^A (21A-'^ + 252^2;^^^ _ 35^2+ 
135£;^A - 24A + 1S5E^ - 252^^ ^ ^S) M'^ + Sr^X^ {21X^ + 3UE^X^- 
95X^ + S60E*X - SiOE^X + lOOA + 360£;* - GSiE^ + 24) + 2rlX^ ■ 
{88E^X^ - 47X^ + ISOE'^X - 260E^X + 25A + 180£;'' - 348E^ - 24) 
+ 2rlX^ {6E^X^ + SOE^X - 53E^X + 23X + 30E'^ 

ArlX^ (A + 1) {E^X - 2E^X - 2) 



59E^ 



11 M + 



(44) 



3KErii 



rl{2M-ruf{3M + ruXf 



135M^ + 27r„ (7A + 32^;^ - 6) + 3r^- 

(35A^ + 396E'^X - 75A + 108^;'' - 144£;2 + 18) + r^A (35A2+ 
612£;2a - 120A + 432^;^ - 594E^ + 72) + r^A^ [UOE'^X - 45A+ 
216^^ - 306£;^ + 36) + 2r^A^ (QE'^X + 24E^ - 35E'^ + 9)M + 



2rlX^ {2E'^X - 3E^X - l) 



(45) 
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-kE 



1431M'^ + 6r„ (251A + 234:E^ - 210) M"* 



9rl (59A2 + mOE^X - 148A + 36£;'' - 66E^ + 30) Af^ + 6r^A (10A2 + 
82£;2A - 79A + biE'^ - 102E'^ + 48) + 2ri\^ {28E^X - 27A + biE'^ 

, 2 



-105i;2 + 52) M + 12r^A^(£;2 _ i)^ 



(46) 



nEf,. 



ri (2M - r„)2 (3Af + r„A)' 



243^/"* + 3r„ (61A + 132^;^ - 64) + 3r^ • 



(12A2 + 92£;2A - 49A + 36^;'* - ^?>E^ + 12) + 2r^A {2AE^\ - 15A+ 
36£;^ - hlE"^ + 14) Af + 6rf,A2 (i;^ - l) i2E'^ - l) (47) 



-kE 



rl (3A// + r„A) 
QE"^ - 5) Af + l2rlX{E'^ - l)^ 



189A/^ + 2r„ (36A + 84£;2 - 77) Af ^ + 6r^ (i;^ - l) (lOA- 



(48) 



While heuristic arguments [35j|36j have been put forward to show that, for radial 
fall in the RW gauge, even metric perturbations belong to the C° continuity class 
at the position of th particle, in 20 21 we have provided an analysis vis d vis the 
jump conditions that the wave function and its (first and second) derivatives have 
to satisfy for guaranteeing the continuity of the perturbations at the position of the 



particle. Therein, we have derived the same jump conditions ( 34 - 38 1 from the inverse 



relations (expressions giving the perturbations as function of the wave function and 
its derivatives) by fulfilment of the continuity conditions (equal coefficients for the 
two Heaviside distributions, and null coefficients for the Dirac distribution and its 
derivative) . 



4. The algorithm 

The integration method considers cells belonging to two groups for cells never crossed 
by the world line, the integrating method may be drawn by previous approaches 
explored by Lousto 



22 and Haas 15 



whereas for cells crossed by a particle, we 
propose a new algorithm. The grid is in the r* , t domain. 

Initial conditions require knowledge of the situation prior to t = 0. At fourth 
order, the wave function may be Taylor-expanded around t — 0. For the boundary 
conditions, simplicity suggests a sufficiently huge grid to avoid unwanted refiections. 



4.I. Empty cells 

Empty cells are those cells which are not crossed by the particle. In this case, the 
cell upper point is obtained by performing an integration of the wave equation over 
the entire surface A of the cell, identified by the nodes a,/3,7,(5. We briefiy recall the 



algorithm used by Haas 15 . Therein, the sole numerical computation to be carried 
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out is represented by the product of the potential term and the wave function V'i!— g. 
It is performed via a double Simpson integral, using points of the past light cone of 
the upper node a, Fig. [T| We set = g{r*,tq) = V{rq)'^{r*tq), Vq = V{rq) and 
^q = \l/(r*,tq), where q is one of the points shown in Fig. flj The increment h is 
defined as h — |Ar* — ^At where Ar* is the spatial step and At is the time step. 




Figure 1: Set of points (circles and crosses) used for the integration of the V'i'= g term in 
the vacuum case. The crosses don't overlap with grid nodes; thus the field g at these points, 
Eqs. ( |50|51 1, is approximated by the field at the nodes on the past light cone of the grid 
node a. 



We have 



gdA^ 



Cell 



5a +5,3 +57 +55 +4 (5/37 + 5q/3 +557 +5a5) + 165(7 +0{h^) , (49) 



where the sum of the intermediate terms between nodes is given by 



507 + 5q/3 + 557 + 9aS = ZK^cr 



^^57*<5 



1 - 



1 rh 



2 \2 

1 /h 



2 V 2 



d'-y 





















1 1 










1 + 2 









0{h^) . 



(50) 



The last intermediate term g^ in Eq. 49 is evaluated using given nodes in the past 
hght cone of a. Fig. [T] 

85/J + 857 + 8.g5 - 4^^^ - 45^2 + .g^j - .g^^ - + + 0{h'^) . (51) 



9a 



16 



For the differential operators, an exact integration simply leads to 



Cell 



{dl, -dl)^{r\t)dA=^-A 



(52) 
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Finally, we get 



1 

1 - - 

4 



1 - 



1 



1 



^^)+16 



+ Vs) + 



1 

16 



5/37 + .9a^i + .957 + 3a5 + 45o- 



(53) 



For cells adjacent to cells crossed by the particle, the requirement of good accuracy 
suggests a different dealing for the computation of ga, since the past light cone of an 
adjacent cell can cross the path of the particle. In such a case, g^ is approximated by 



non-centred spatial finite difference expressions 15 



4-2. Cells crossed by the world line 

For a given cell, our aim remains the determination of the wave function value at the 
upper node, now rebaptised ao. As in the previous section, we consider (fifteen) points 
both located in the past light cone of the ao point and lying around a chosen point on 
the discontinuity r^it), with the intent of determining ^>ag by their linear combination. 
The non-regularity of the wave function due to the discontinuity, obviously entails a 
different value according to whether the discontinuity is approached from below (^~, 
left of the trajectory. Figs. [2]-[4| or above {^~^, right of the trajectory. Figs. [2]-|4]) 
the particle in radial fall. The same stands for the wave function derivatives. The 
addition of the jump condition to the value of the e.g. (^^) wave function (or 
derivative of) allows to equate this sum to the value '5+ (^~) of the wave function (or 
derivative of). This straightforward property turns being helpful for the achievement 
of the just mentioned linear combination of fifteen points. Incidentally, other linear 
combinations may be envisaged, though combinations of points located solely on one 
side of the discontinuity are to be avoided. 

With reference to Figs. [2]-|4j there are three different cases depending upon how 
the trajectory of the particle crosses the cell wherein ao lies. These three cases are 
further subdivided into three sub-cases, for a total of nine. In the following, we label 
by R the points on the right of the [aoae] line and by L the points on the left. Dealing 
with radial fall, and thereby with a 2D code, the up and down labels might be proper; 
nevertheless, we stick to right and left labels, given the orientation of the r* axis in 
the Figs. [2] -[4] For the first group of three, the trajectory crosses the [a2/3f^] and 
[ao/Sf] lines, Fig. [Sj for the second group, the [a2/3f ] and [ao/Sf] lines, Fig. [sj finally 
for the third group, the [a2/3f^] and [ao/3f ] lines. Fig. [i] 

We start considering the sub-case (la) shown by Fig. [5] for which the trajectory 
crosses the line [aoa2] at the point b. For compactness of the presentation of the 
final results, while we still adopt the same notation for the jump conditions, namely 
[^]q for the difference ('5+ — ^'~)^ (j ), for the jump derivatives instead, we rely 

henceforth on the notation [a,".^^*], = (a;?.5t'"«'+ - i9;-i9r*")r„=r„(t,)' where tg is 
the coordinate time at the point q = a,b. We also define the lapse £{, = ta^ —ti,. 

We recall that our aim is the determination of the value of '^^o^ knowing: i) ei,, ii) 
the jump (analytical) conditions on '3/ and its derivatives at the point b; iii) the values 
of \1/ on a set of fifteen points {a, 7, fi, v} at the left and right sides of the world 
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line. A Taylor series is applied at each point around b up to fourth order, thereby 
obtaining 

K.^'^^{h + et,rl)^J2'i^,^t'ft + 0{el) , (54) 
^a.-^"fa-(^fe-^b),^.*)-E(-l) ^, dl'^t+0{h') , (55) 

^^^n^,^^^{tt-{jh-et),r;±h)^ E (_1)'"(±1)" 1:^ J^OnQrn^±+o{h^) , 

n+m<4 

(56) 

f± =*±(4-(fc/i-e.),r:±2/.)= (-1)'"(±1)" ^^T ^^^~f"^ d-,dr^i+0 {h') , 

ifc 71. TTl. 

n+m<4 

(57) 

*^K..=*^(i.-(3/i-e.),r,*±3/i)= V {-ir{±ir^^i^hz^ongn.^±+o{h') , 
Ms ■'^ — ^ n! m! 

n+m<4 

(58) 

4'\,.-**(i6-(4/i-e,),r,*±4/i)= V (_l)™(±l)»il^il^iz£^a« a^vpi+O (/,5) , 

n+7n<4 

(59) 

for the indexes running as i = 2, 4, 6, j = 1, 3 and A: = 2, 4 and concerning the a, (3 and 
7 nodes, respectively. Our notation implies that the subscript R, L stands for R when 
the superscript ± corresponds to +, whereas R,L stands for L when ± corresponds 
to — . With reference to Eq. [54j we get 

4 4 

- E C"^"^^ + ^ (^') = E ^» (^"*^ + 1^"*];.) + ^ e^') = 

4 

= co*^ + c,dt% + c^a^M/- + C39f + c,dt% + E cn [ar*], + O {h') (60) 

n=0 

For an accuracy at fourth order, all quantities 0{h^) are disregarded. The 



sum S = co*b + ci^t^'f, + C2df'iff, + c^df^i/^ + c^df'^f^ , Eq. 60 is composed 
by numerical derivatives of lower order than 0{h^), and therefore they can't be 
neglected. However, the computation of high order derivatives is often accompanied 
by numerical noise. Therefore, we replace this sum by a combination of wave 
function values in the ao light cone. This is attained in two steps. The 
former involves taking fifteen wave function values on the two sides of the 

trajectory, that is , ^-"^ , ^'+«, , *+„, , , }^ Fig. jij The 

latter employs the jump conditions to relate the fifteen mentioned points with 
I*- For the former step, we define the 
sum S 
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+ M^^, + Mf^+ + + M^^+n ■ (61) 

^■3 f^3 ^4 ^4 

where {A„ Bf , B^,gj^ ,gj^ , ,M§,Afi' ,Aff} are constants. 

We observe that the 5 sum entails only wave function values at the left of the 
6 point on the trajectory. The jump conditions are once more exploited to relate 
the two domains r* < rl^{t) and r* > r*(t). This specifically concerns six points 
7^, z/|*^}. For instance, at the point, we can write 




Figure 2: The three sub-cases for which the particle enters through the [02/31*] side and 



leaves through the [ao/3i ] side. The elimination of the ^'^ derivatives demands, Eq. 60 
the utilisation of fifteen points, represented by circles, in the light cone of qq- Numerical 
efficiency suggests that the points are taken at both left and right sides of the r'^{t) 
trajectory. In the three cases, the particle crosses the line [0002] at the point b. The 
background distinguishes two zones: one where 3'(r* <r*(t),t) = "^'{r* ,t), the other where 
^'(r* > r^{t),t) = '^^{r*,t), the path r^{t) representing the separation between the two 
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E (-i)"^ J {d-,dr^- + [d-,dr^,)+o{h'') 

n+m<4 

= ^"«+ E (-ir^^^^^^[5;iar*], , (62) 

n\ ml 

n~\-rrL<4 

where 

%f= E (-ir^^'^;^(a"5r*^)+^(^') • (63) 

n+m<4 





Figure 3: The three sub-cases for which the particle enters through the [a2/3f ] side and leaves 
through the [oo/Sf ] side. The elimination of the 5'^ derivatives demands the utilisation of 
fifteen points, represented by circles, in the light cone of ao. Numerical efficiency suggests 
that the points are taken at both left and right sides of the r^(t) trajectory. In the three 
cases, the particle crosses the line /^i*] at the point a. The background distinguishes two 
zones: one where 4'(r* <r*(f),f) = ^'~(r*,f), the other where ^{r* >rl^{t),t) — 'I'"'"(r*, t), 
the path r^{t) representing the separation between the two zones. 
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(3c) 



Figure 4: The three sub-cases for which the particle enters through the [02/31*] side and leaves 
through the [qo/3i*] side. The elimination of the ^I*" derivatives demands the utilisation of 
fifteen points, represented by circles, in the light cone of ao. Numerical efflciency suggests 
that the points are taken at both left and right sides of the r*(t) trajectory. In the three 
cases, the particle crosses the line /Sf*] at the point a. The background distinguishes two 
zones: one where '^{r* <r^{t),t) = '^~{r*,t), the other where '^{r* >r^{t),t) — 4'^(r*,t), 
the path r* {t) representing the separation between the two zones. 

By application of the same transformation to the quantities , 'i'~^R, Eq. 

Tfc ^3 "3 

161] becomes 



i j k 

+ A^f * . + Mi^„ + AAf + A(f , (64) 

A^3 M3 ^4 ^4 

where <i>^"™^ is an analytic function, composed by the jump conditions at the b point, 



weighted by coefficients issued by Eq. 62 or similar equations 



Having only 5" terms on the right hand side of Eq. |64j we can finally search 
the coefficients {Ai, Bf , ,0^ ,Q^,M3 , M§,Af4 ,J\fi^} that satisfy the equation 
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+ M^^r + Mi^„ + N^%, + . (65) 

^3 ^3 I/^ I/^ 

Using the notation of Eqs. |62] [63j and by injection of Eqs. [55]|59l a Taylor 



expansion of fourth order at the b point is apphed to the right-hand side of Eq. 65 
The system can be cast in a matrix form 

T • P = C , (66) 

where P is the unknown 15-vector formed by the coefficients {Ai, ^Bf, Gk^Gki -M3 , 

P^{A2,A4,A6,B^,Bi,B^,B^,gi,gtg^,G^,M^,Mi,Ml,Mfy , (67) 

and C is given by the 15-vector 

C = (co, ci, C2, C3, C4, 0, • • • , 0)* , (68) 

while T is the (15 x 15) matrix constructed from the Taylor coefficients in Eqs. 55]|59 
(see appendix). By inversion of T, we get P and specifically 



A2 = 


-27 


A4 = 


-9 1 
5 ' " 5 ' 


B^ = 


6f = 


12 

y ' 




G^ = 


5? = 


-9 


rL rR 


■M3 = 


= Mi 


2 

" 5 ' 


= AAf = 



The following equivalences path the last stretch of the way 



4 



= 5 - + ^ c„ [ar*], = s + $1,^.) , (69) 



n=0 



and explicitly, we get 



27 ^_ 9 ^_ 1 _ 12 /_ ^ , \ 18 



where ^'Jh = ^^^^ for sub-case (la), and ^^^^ = for sub-cases (lb,lc); = ^^p. 

"3 P3 ^3 ^3 T4 ')'4 

for sub-cases (la, lb), and '^^r — for sub-case (Ic); and is an analytic 

function, that for the (la) sub-case, it takes the value 
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lb [^**]^ 

5efc4_56/ieb3^216/i2eb2_ 384/j3g^_|.240/i4 13/1.^ . 

i^t - — [^--^Ift 

(71) 

The quantity varies according to the different sub-cases: for the case (lb) 
of Fig. [2] the point 13^, whereas for the case (Ic) the points and 74" are in the 
r* < r* domain. Therefore 

£b'^12/teb^ + 54/i^ efc-66/t3 

10 L * Jb 

- 16 /i efo^ + 108 _ 254 246 /i" .^4 . 
+ i5 [^**]& 



(72) 



+ '-^^^^ [dl-^d.^,~ 'JIS^ [ea.^], 

^— ^ ^ [dr-dM, ■ (73) 

We thus have obtained, without direct integration of the singular source and the 
potential term, the value of the upper node. The equations shows three types of 
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terms: the preceding node values of the same ceU, the jump conditions which are fuhy 
analytical quantities, and the wave function values at adjacent cells. Incidentally, 
at first order f21l, the latter type of terms disappears and a simpler expression is 
obtained. 

Similar relations are found for the other two remaining cases. For case 2, Fig. [3j 
we obtain (having defined the shift = tpR—r*) 

5 5 5 5 V '3f + ^/3f J + 5 V + */33« 

9 - ^ ^ 3 + \ 2 / , \ (2) 



as ae 
(2) 

'■:(ta) 



where = for sub-cases (2a,2b), and = ^'^^ for sub-case (2c); ^"^^ = 5-+ 
for sub-case (2a), and ^'^ 
takes the following value 



for sub-case (2a), and ^'^ = for sub-cases (2b, 2c). For the (2a) sub-case, <1'*'^'' 



, 4 (5ea - 8/.) ^^^^^j^ ^ 2(6.-fe) (5 6, -11/.) ^^^^^^ 



5 ^ ' 'a ■ 5 

^ ^5 t^'-^l 

(e, -/i) (5e,3-27/ie,2+39/i2^^_5;j3) 



30 



- [d..d,^l - M^g-/^) (11^.-23/.) 

+ ^^^^^ [^5.^*], - [a.a3*]^ . (75) 

For the same preceding reason, the sub-cases (2b, 2c) differ as the points a4 and 
ag are or aren't in the r* > r* domain. Therefore, we have 



(26) 42,^, 27/i,^^, 69/i2 2,1 13/i3.3 . 231 /i'^ 



3 (7ea - 11 h) ^^^^ , 3 (e, - h ) (76,- 15fe) ^^^2 

10 

(e,-/t) (7ea3-37/i£a^+53;i^£, -7fe3) 4 
40 ^ 

3h^ (236. - 31/.) ^^^^^2^^^ _ Shje^^hf (36.-7/.) ^^3^^^^^^ 



10 L ' t Ja 10 



h^ (eg-h) (1576. - 109 /.) . 2 ^2^] _ fe' (65 6.- 69 /i) ^3 
20 L * J" 10 
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(76) 



^ehHe hf ^h^ ie.^h)il9e.^7h) 

,77, 

Finally for case 3, Fig. |4j we have 

27 ^_ 1_ 12 /_ ,\ 18 /_ 



where ^'^h = for sub-case (3a), and — for sub-cases (3b, 3c); = 

Ps P3 P3 P3 74 74 

for sub-cases (3a, 3b), and ^I^^ji = ^f^j, for sub-case (3c); and ^ takes the values 



1(1) - + — [a,*]^ - — [d^^]^ + — - — 



2 (5e„3-6fe£,2-3/i^ea + fe3) 3 ^ 

15 L Ja 

(ea-/i) (5ea'-3/ie,2_9/i2e^ _5/i3) 



30 L'^-'^Ja 



2/. (ll£„-5fe) ^^^^ ^ ^ (6.-fe) (ll6, + fe) 



2/.^ (ll.„-5/.) [g^,,.^]^^Me.-M- (ll.. + 7M [,3^,^^]^ 



5 L ' I Ja ■ 15 



+ + (35. -29/.) ^^^^^3^^^ ^ ^^^^ 



(36) 2. . 14/1,^-,, 14/i\o2Tn 37/13 3 ii;j4 
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15 t^-*]'^ 

ea^ -I6h - 12 /l^ 4, g ^3 _^ /l^ 

60 

5 5 
2h^{7ea + 5h) h {7ea^ + l5hea^ + 3h^ea~7h^) 

5 L * J a L !•* t J Q 

-^^ (7^.^ + 10^^. + ^^) (37 .. + 29/.) ^ 

+ (/. - _ 5eJ - 10 he.-h^ 

30 

5 - 20 /i - 6 /i^ e„2 + 28 h^ + 23 /i" 
120 

The jump conditions in the tortoise r* relate to those previously computed in the r 
variable (the relations for mixed derivatives {r*,t) are easily inferred) 

= /,■„ , (82) 

[^,r'r'] = frJi .r] + /'„ [*,r.] , (83) 

- /r„ + + ^flfr^ ,rr] + .rrr] , (84) 

+ 6/3 [f^rrr] + fr^ [^-.rrrr] ■ (85) 

5. Numerical implementation 

Waveforms at infinity and at the particle position at first order are to be found 



in 21 , as well as comparisons with other methods. Herein we are concerned on 
the numerical improvement. To this end, we have considered a distant observer, 
located at r* = 400(2A/). The observer is reached by a pulse produced by a Gaussian, 
time-symmetric perturbation 



*(r*,t)t=o = exp[-(r*-rS)2] , (86) 
dt^{r\t)t=Q = . (87) 
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Fig. [5j obtained for r„o — 5(2M), shows the waveform produced in the 
homogeneous case. The convergence rate is computed as (e'-"-'(0 is the unknown 
error function of order w 1) 




400 420 440 460 480 



t/2M 



Figure 5: The waveform, ruo = 5(2M), of a Gaussian, time-symmetric initial pulse. The 
observer is located at r* = 400(2M). 




Fig. [6] obtained for r„o — 5(2M), shows the fourth and second order convergence 
rates (we remind that the first order code ^l] includes empty cells dealt at second 
order) . 
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6. Conclusions 



We have presented a fourth order novel integration method in time domain for the 
Zerilli wave equation. We have focused our attention to the even perturbations 
produced by a particle plunging in a non-rotating black hole. For cells crossed by 
the particle world line, the forward time wave function value at the upper node of 
the {t, r*) grid cell is obtained by the combination of the preceding node values of the 
same cell, analytic expressions related to the jump conditions, and the values of the 
wave function at adjacent cells. In this manner, the numerical integration does not 
deal directly nor with the source term and the associated singularities, nor with the 
potential term. In short, the direct integration of the wave equation is avoided. For 



other cells, we refer instead to already published approaches 15 



The scheme has also been applied to circular and eccentric orbits and it will be 
object of a forthcoming publication. 
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Appendix 

Through Eq. |60j we have determined the value of 5* at the upper node of the cell 
as function of the analytic jump conditions and of the time derivatives of the wave 
function up to fourth order. The derivatives are evaluated at the point b and weighted 
by five coefficients cq, ci, C2, C3 and C4. Afterwards, the derivatives are converted into 
a linear combination of the wave function values taken on points at the left and right 
sides of the trajectory. Indeed, Eq. |65] represents such a system of linear equations. 
By injection of Eqs. [56p9| into Eq. |65j we get 

+ 

+ 
+ 

where Tp"'"^^ represent the Taylor series coefficients at p in the neighbourhood of b 
and the indexes correspond to n*'' space and to*'' time derivatives. The wave function 
at p is thus given by 

= J2 T'^"'™^<9;'.ar*^ + ■ (90) 



+ 
+ 

= 0292*- ,(89) 
+ 

+ 
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An exam ple shows the procedure which is appHcable to all cases. We pick the node 02, 



Eq. 



55 



where Ta 



(-1) 



and remind that T, 



(0,0) 



1 V p. By grouping 



the derivatives, we get 



(0,0) 



(0,0) 



.••+A/-4Vr)vI;- 



+ 
+ 

+ 



(91) 



By identification, we obtain a linear system, that is cast in the form 

1 \ 



/ 1 • 

7^(0,1) 


1 

7^(0,1) 

■ hi ■ 


1 

7^(0,1) 


1 

y(0,l) 


7^(0,4) 


7^(0,4) 

■ hi ■ 


7^(0,4) 


y(0,4) 


7.(1,0) 


■ hi ■ 


7.(1,0) 


2.(1,0) 




7^(1,3) 


7.(1,3) 


2.(1,3) 


T 



.(0,1) 



.(0,4) 



7.(1,0) 



T 



(1,3) 



/ ^2 \ 



/ '^o \ 

ci 

C2 
C3 
C4 





V j 



where the upper indexes (n,m) cover all combinations such that n + m < 4. Finally, 
by inversion of the T matrix, the unknown terms of the P vector are identified. 
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